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1  Introduction 

1.1  Motivation  for  Research 

Naval  architecture  in  the  United  States  Navy  has  come  full  circle.  With  the 
selection  of  a  tumblehome  hull  for  the  DDG-1000  project,  the  United  States  Navy 
has  come  back  to  a  design  concept  that  was  used  on  the  oldest  commissioned 
warship  in  their  fleet,  the  USS  Constitution.  Tumblehome,  the  inward  slope  of 
the  upper  part  of  the  hull  of  a  ship,  has  come  back  in  style  due  to  the  stealthy 
benefits  it  imparts  upon  its  ship.  However,  with  the  use  of  this  ‘new’  concept 
come  interesting  challenges  in  naval  architecture. 

The  latest  designs  for  the  hulls  of  United  States  naval  warships  have 
sparked  a  controversy  over  their  stability.  In  the  April  2007  edition  of  the  Navy 
Times,  experts  are  quoted  stating  that  the  hull  may  be  unstable  in  roll  motion. [1] 
Programs  that  use  linear  sea  keeping  approximations  predict  good  stability  for 
tumblehome  hulls.  However,  non-linear  effects  in  the  coupling  of  pitch  and 
heave  with  roll  motions  can  degrade  the  stability  of  a  tumblehome  ship. 
Programs  that  calculate  the  non-linear  effects  of  this  nature  may  require 
significant  computing  power  making  them  costly  and  time  consuming.  A  simpler 
approach  could  reduce  computing  requirements  and  allow  naval  architects  to 
quickly  estimate  the  stability  of  any  new  hull. 

1 .2  Objective  and  Outline  of  Thesis 

The  purpose  of  this  thesis  is  to  explore  the  non-linear  effects  that  pitch, 
heave,  and  roll  exert  on  a  ship’s  restoring  moment.  The  linear  seakeeping 
equations  will  be  modified  to  allow  the  restoration  force  to  be  calculated  as  a 
function  of  roll,  pitch,  and  heave.  Time  simulations  of  pitch  and  heave  responses 
derived  from  linear  seakeeping  theory  are  used  to  provide  a  basis  for  numerical 
solution  of  the  non-linear  seakeeping  equation  for  roll.  Although  the  research 
presented  is  for  one  specific  hull  form  and  operating  condition,  the  methods  used 
are  intended  to  applicable  to  any  ship  in  any  seakeeping  situation.  Chapter  Two 
discusses  the  characteristics  of  the  hull  used  in  this  research,  including  the 
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operational  scenario  that  is  used  for  conducting  the  research.  Chapter  Three 
explains  linear  seakeeping  theory,  as  well  as  the  approximations  used  by  the 
seakeeping  program  Maxsurf.  Chapter  Four  discusses  the  method  used  to 
determine  the  excitation  forces  used  in  the  non-linear  simulation.  Chapter  Five 
describes  the  computer  programs  used  to  simulate  the  roll  motion  of  the  ship. 
Chapter  Six  presents  a  comparison  of  the  linear  and  non-linear  results,  and 
Chapter  Seven  presents  recommendations  for  future  work  and  conclusions. 
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2  The  ONR  Tumblehome  Hull 


2.1  Hull  Description 

The  hull  form  used  for  this  thesis  is  the  tumblehome  hull  used  by  the 
Office  of  Naval  Research.  This  ship  has  a  length  of  154  meters  and  a  beam  of 
18.5  meters.  The  displacement  at  a  draft  of  5.5  meters  is  approximately  8800 
metric  tons. [2]  The  hull  has  a  vertical  center  of  gravity  assumed  to  be  7.5  meters 
above  baseline,  and  the  transverse  metacentric  height  is  1.47  meters.  The  bow 
rakes  aft  instead  of  forward  necessitated  by  the  tumblehome  sides.  The  profile 
and  body  plan  of  the  ship  is  shown  in  Figure  2-1  and  Figure  2-2: 


Figure  2-1 :  Profile  View  of  the  ONR  T umblehome  Hull. 
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Figure  2-2:  Body  Plan  of  ONR  Tumblehome  Hull 

The  tumblehome  hull  has  military  advantages  that  make  it  attractive  for 
use  in  surface  combatants.  The  chief  advantage  comes  from  the  fact  that  the 
sides  of  the  hull  are  angled  away  from  the  waterline.  This  will  tend  to  reflect 
radar  energy  that  is  directed  towards  the  ship  from  another  up  into  the 
atmosphere,  and  not  back  towards  the  radar  receiver.  This  scattering  of  radar 
energy  significantly  reduces  the  radar  cross  section  of  the  ship,  making  the  ship 
harder  to  detect. 

From  a  naval  architecture  standpoint,  tumblehome  hulls  have  some  less 
than  ideal  properties.  The  decreasing  beam  above  the  waterline  limits  the 
deckhouse  area,  and  reduces  the  damaged  stability  of  the  hull  form  compared  to 
wall  sided  or  flared  hulls.  In  the  area  of  stability,  the  hull  suffers  compared  to 
traditional  hull  forms.  As  the  ship  lists,  the  tumblehome  causes  a  relative 
reduction  in  the  waterplane  area  of  the  ship,  lowering  the  righting  arm.  The 
righting  arm  for  zero  pitch  is  shown  in  Figure  2-3,  compared  to  ships  with  similar 
shapes,  but  with  wall  sides  and  flare  instead  of  tumblehome. 
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Medium  GM  (1.5m)  Righting  Arm  Curves 


Figure  2-3:  Righting  Arm  Curves  for  Different  Hull  Shapes 

The  wave  piercing  bow  has  its  own  military  advantages.  Because  of  the 
reduced  buoyancy  of  the  bow,  the  ship  tends  to  pitch  less  as  it  rides  through 
waves,  giving  the  ship  a  more  stable  platform  for  its  weapons  systems. 
Additionally,  the  design  may  reduce  the  waves  generated  by  the  ship,  lowering 
resistance  of  the  hull.  These  advantages  come  at  a  price  in  the  stability  analysis. 
The  ship  loses  waterplane  area  as  it  pitches  down  by  the  bow,  resulting  in  a 
smaller  transverse  righting  arm  as  the  bow  sinks  into  the  water.  In  most 
operating  cases,  this  effect  can  be  neglected.  In  heavier  seas  where  pitch  is 
significant,  the  loss  of  righting  arm  can  cause  severe  rolling  transients. 

2.2  Computer  Modeling  in  Maxsurf  and  Matlab 

The  ONR  tumblehome  hull  was  provided  by  the  Carderock  Division  of  the 
Naval  Surface  Warfare  Center.  In  order  to  conduct  this  research,  it  was  modeled 
in  both  Maxsurf  and  Matlab.  Offsets  were  generated  in  Maxsurf  using  a 
computer  aided  drawing  of  the  hull.  The  offsets  were  used  in  three  different  parts 
of  this  research. 

First,  Maxsurf’s  Sea  Keeper  module  used  the  offsets  to  determine  a  pitch, 
heave,  and  roll  response  using  the  program’s  linear  strip  theory  calculations. 
These  outputs  were  in  the  form  of  linear  response  amplitude  operators  (RAOs) 
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for  both  magnitude  and  phase  relative  to  the  incident  waves.  As  explained  in 
Chapter  Five,  the  RAOs  were  used  to  generate  a  time  simulation  of  the  ship’s 
pitch,  heave,  and  roll.  The  pitch  and  heave  time  simulations  were  used  to 
determine  the  non-linear  roll  righting  moment.  The  roll  response  was  used  to 
compare  against  the  results  of  the  non-linear  simulation.  The  offsets  were  also 
used  to  model  the  ship  in  Matlab  to  determine  excitation  forces.  The  offsets  were 
used  as  a  set  of  nodes  to  compute  the  Froude-Krylov  and  diffraction  excitation 
moments.  Finally,  the  hull  offsets  were  used  in  Maxsurf’s  Hydromax  module  to 
determine  the  ship’s  righting  moment  in  roll  as  a  function  of  pitch,  heave,  and  roll. 
This  will  be  explained  in  Section  2.3. 

Mass  distribution  for  the  ship  was  chosen  to  approximate  a  ship  with  a 
relatively  high  center  of  gravity  to  allow  for  high  weight,  such  as  weapons 
systems  on  a  naval  vessel.  The  center  of  gravity  was  chosen  to  be  on  the 
centerline  two  meters  above  the  still  waterline.  This  was  arbitrarily  chosen  so 
that  the  metacentric  height  would  match  the  height  shown  in  Figure  2-3.  The 
longitudinal  position  was  determined  as  directly  above  the  center  of  buoyancy  to 
give  the  ship  no  trim  in  pitch.  In  order  to  match  Maxsurf  approximations,  an 
estimate  of  40  percent  of  the  ship’s  beam  was  used  to  approximate  the  ship’s 
radius  of  gyration.  Therefore  the  roll  moment  of  inertia  is  given  by  the  formula: 

/,4=(0.45)VV  (2.1) 

This  mass  distribution  may  be  used  to  estimate  the  undamped  natural 
period  of  the  hull  [5]  through  the  use  of  the  equation 


T  =  2n  *^^44  +  ^4)  (2.2) 

V  p^gGM, 

The  added  mass  in  roll,  A44,  was  assumed  to  be  0.3  times  the  mass 
moment  of  inertia,  as  discussed  in  Section  3.4.  For  the  values  given  in  this 
section,  the  natural  period  is  approximately  13.9  seconds.  This  will  be  used  to 
compare  against  results  in  Chapter  Six. 
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2.3  Non-Linearities  in  Righting  Moment 


One  of  the  unique  aspects  of  a  tumblehome  hull  is  its  roll  restoring 
moment  behavior  with  its  movement  in  the  three  principal  degrees  of  freedom. 
Because  of  the  wave  piercing  bow,  its  buoyancy  forward  is  much  smaller  than  its 
buoyancy  aft.  This  asymmetry  leads  to  a  non-linear  behavior  in  restoring 
moment.  As  mentioned  in  the  previous  section,  the  offsets  were  used  to 
determine  the  hydrostatic  righting  moment  as  a  function  of  roll  angle,  pitch  angle, 
and  draft. 

Righting  moment  varies  significantly  with  both  pitch  and  draft.  In  the  linear 
simulations  determined  from  Sea  Keeper,  pitch  varied  between  15  meters  to  -15 
meters  by  the  stern,  and  draft  at  amidships  varied  between  3  meters  and  7.5 
meters.  This  behavior  is  shown  graphically  in  Figure  2-4  and  Figure  2-5.  Wave 
elevation  and  slope  effects  are  not  included  in  the  restoring  moment 
determination.  These  effects  were  left  for  objects  of  further  study. 

3.5 

3 

?  2.5 

z 

c 

oi 

o  2 

2 


0.5 

20  15  10  5  0  -5  -10  -15  -20 

Trim  (+  is  by  stern)  (m) 

Figure  2-4:  Righting  Moment  vs.  Trim  for  Roii=10,  Draft=5.5m 
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Figure  2-5:  Righting  Moment  vs.  Draft,  Roii=10,  No  pitch 

Trim  effects  were  most  drastic  near  zero  trim.  As  the  ship  pitches  by  the 
bow,  the  righting  moment  decreases  as  expected.  This  leads  to  an  obvious 
design  technique  for  improved  sea  keeping  of  trimming  the  vessel  by  the  stern. 
However,  this  method  of  improving  the  righting  moment  reduces  the  advantages 
of  having  a  wave  piercing  bow  in  the  first  place.  The  draft  effects  are  not  as 
significant  as  the  effects  for  pitch,  but  a  low  draft  coupled  with  a  significant  pitch 
by  the  bow  could  lead  to  very  small  righting  moments.  In  very  extreme  cases, 
the  moment  could  turn  to  an  overturning  moment.  Obviously,  this  situation 
should  be  avoided  at  all  costs. 

2.4  Seakeeping  Scenario 

The  seakeeping  scenario  chosen  is  well  outside  a  normal  operating 
envelope  for  naval  vessels.  The  conditions  were  chosen  to  excite  significant 
motion  in  any  ship.  The  ship  was  set  to  a  speed  of  25  knots,  traveling  in  Sea 
State  Eight  with  seas  30  degrees  abaft  of  the  beam.  The  significant  wave  height 
for  Sea  State  Eight  is  11.5  meters,  with  a  modal  period  of  0.37  radians  per 
second.  The  seas  are  assumed  to  be  fully  developed,  long-crested,  waves  from 
a  distant  storm.  A  sketch  of  this  scenario  is  shown  in  the  figure  below: 
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A  consistent  coordinate  system  was  chosen  for  all  parts  of  the  study.  The 
coordinate  system  is  the  Cartesian  coordinate  system  normally  used  for  sea 
keeping  problems.  The  x-axis  is  directed  along  the  longitudinal  center  of  the 
ship,  and  the  z-axis  is  pointed  so  that  positive  z  is  up.  Therefore,  the  y-axis  is 
directed  out  the  port  side  of  the  ship.  The  six  degrees  of  freedom  are  therefore 
surge  in  the  x-direction;  sway  in  the  y;  heave  in  the  z;  roll  about  the  x-axis;  pitch 
about  the  y;  and  finally  yaw  about  the  z.  The  origin  was  selected  as  the 
centerline  at  amidships,  at  a  height  equal  to  that  of  the  still  waterline. 

For  the  seakeeping  part  of  the  problem,  roll  was  assumed  to  be  decoupled 
from  sway  and  yaw.  While  this  is  not  the  actual  case  as  discussed  in  Section 
3.3,  it  provides  a  simplification  to  the  problem  to  get  a  basis  for  comparison  with 
another  seakeeping  program. 

3  Linear  Seakeeping  Theory 

3.1  Linear  Plane  Progressive  Waves 

Inside  the  domain,  random  seas  were  simulated  by  unidirectional  plane 
waves  incident  upon  the  hull  at  an  angle  p  from  the  stern  of  the  ship.  The  plane 
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waves  were  assumed  to  be  a  superposition  of  sinusoidal  waves  at  different 
frequencies.  To  treat  the  waves  using  potential  theory,  the  following 
assumptions  were  made.  First,  the  unsteady  viscous  forces  on  the  hull  are 
neglected.  Next,  the  density  of  the  seawater  is  constant  and  all  flows  were 
incompressible.  Finally  flow  was  approximated  as  being  irrotational. 

All  flow  potential  functions  must  satisfy  conservation  of  mass.  With  the 
assumption  of  incompressible  flow,  the  conservation  of  mass  states  that  the 
divergence  of  the  velocity  field  must  be  zero.  In  mathematical  terms,  this  means 
that  the  governing  equation  for  all  potentials  must  be: 


du' 

dx'  dz'  dx'^  dz^ 


(3.1) 


In  equation(3.1),  the  primed  quantities  are  with  respect  to  a  coordinate 
system  with  the  waves  propagating  in  the  x’  direction.  The  two  boundary 
conditions  for  these  come  from  the  kinematics  of  the  waves  and  the  dynamic 
pressure  of  the  wave’s  surface.  The  kinematic  boundary  condition  is  that  along 
the  wave  surface,  z=^  where  (x,y,t)  is  the  wave  elevation  as  a  function  of 
position  and  time.  By  defining  a  function  F=z-  ,  then  the  kinematic  boundary 
condition  may  be  written  as: 


—  =  — +  (F.V)F  =  0  onz  =  ^  (3.2) 

Dtdt 


The  dynamic  free  surface  condition  comes  from  the  fact  that  the 
pressure  is  atmospheric  along  the  surface  of  the  wave.  From  Bernoulli’s 
equation  for  potential  flow, 

5®  1  p 

—  +  -VO.VO  +  gz  =  — ^  onz=^  (3.3) 

dt  2  p 


In  most  case,  the  wave  slopes  encountered  are  small,  so  these 
equations  may  be  treated  by  linearizing  the  boundary  conditions.  The  linearized 
boundary  conditions  become: 


+  g 


&  ^  =  - 
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1  ao 

g  dt 


on  z=0 


(3.4) 


The  solution  to  this  boundary  value  problem  in  deep  water  with  x 
directed  along  the  direction  of  the  propagation  of  the  waves  in  a  stationary 
coordinate  system  is: 


®  =  (3.5) 

with  the  dispersion  relationship 

o/=kg  (3-6) 

Because  the  waves  are  traveling  at  an  angle  with  respect  to  the 
longitudinal  axis,  the  waves  must  undergo  a  coordinate  transformation.  It  is 
convenient  to  work  in  complex  space  without  the  real  notation,  so  the  complex 
wave  potential  incident  upon  the  ship’s  hull  for  a  given  frequency  is: 

-ik(xcos/^-ysm  J3)+ia)t  kz 

= — e  e  w-/) 


3.2  Wave  Spectra 

Because  real  seas  are  never  monochromatic,  a  statistical  approach  is 
needed  to  simulate  random  seas.  The  approach  used  in  this  thesis  assumes  that 
the  waves  experienced  are  unidirectional,  as  if  they  all  come  from  a  distant 
storm.  It  also  assumes  that  the  surface  elevation  is  a  zero  mean,  ergodic 
collection  with  a  Gaussian  distribution.  The  energy  of  the  storm  may  be 
contained  in  a  spectrum  of  frequencies. 

A  common  spectrum  used  is  the  Bretschneider  Spectrum.  This  is  the 
standard  spectrum  recognized  by  the  International  Towing  Tank  Convention,  and 
will  be  used  throughout  this  paper.  The  Bretschneider  Spectrum  has  the  form: 


S,{G>)  = 


1.25  0)1 


CO 


H,,e 


.\25*{-^y 


(3.8) 
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This  spectrum  requires  two  parameters,  H1/3,  or  significant  wave  height, 
and  modal  frequency,  cOm-  The  significant  wave  height  is  the  average  of  the  one 
third  highest  waves  in  the  seas,  and  the  modal  frequency  is  the  single  frequency 
that  is  most  likely  in  the  distribution.  For  fully  developed  storms,  the  modal 
frequency  is  related  to  the  significant  wave  height  by  the  equation: 


(3.9) 


With  this  relationship,  the  spectrum  is  called  the  Pierson-Moskowitz  spectrum.  In 
Sea  State  Eight,  these  values  are  Hi/3=1 1 .5m,  and  a)m=0.37  rad/s.  A  graph  of 
the  spectrum  is  shown  below. 


Figure  3-1:  Fully  Developed  Bretschneider  (Pierson-Moskowitz)  Spectrum  for  Sea  State 

Eight 
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When  scaled  by  pg,  the  area  under  the  spectrum  curve  represents 
the  energy  per  unit  area  of  the  sea  state.  Because  of  the  assumption  that  the 
sea  state  is  described  by  a  collection  of  waves  at  distinct  frequencies,  this  energy 
relationship  can  be  used  to  simulate  a  random  sea  state  in  a  discrete  form.  For  a 
random  sea,  the  energy  density  of  the  system  may  be  expressed  as 

#  freqs  a2  20 

Ps\s,{(»)d(o  (3.10) 

i=\  ^  0 


By  discretizing  the  integral  in  the  above  relationship,  the  wave 
amplitudes  for  discrete  frequencies,  separated  by  A®,  can  be  calculated  as: 

(3.11) 

Adding  in  a  random  phase  angle,  random  seas  may  be  simulated 
on  a  computer  rather  quickly  using  only  the  formula  for  the  wave  height 
spectrum.  The  surface  elevation  on  a  fixed  (x,y)  vertical  line  ^(r)  is  simply 

#freq 

^ (0  =  X  4  cos(®,r  +  )  (3.12) 

i=\ 

An  example  of  the  record  of  a  random  sea  in  Sea  State  Eight  is  shown  in 
figure  3-2.  The  sea  state  was  generated  using  500  discrete  frequencies  equally 
spaced  from  0  to  2.5  radian/second,  calculated  every  .01  seconds. 
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Figure  3-2:  Simulated  Wave  Elevation  Record  for  Sea  State  Eight 


3.3  Equations  of  Motion  &  Linear  Superposition 

Most  seakeeping  analysis  starts  with  the  assumption  that  the  motions  of  a 
ship  are  linear  in  nature  [3].  This  a  reasonable  due  to  the  small  wave  slopes 
experienced  during  most  normal  operating  scenarios.  The  linearized  seakeeping 
equations  of  motion  for  the  six  degrees  of  freedom  may  be  written  in  short  as 


k=l 


+  Ajk)^k  +  +  C.,?7]  = 


j=l-6 


(3.13) 


The  matrices  M,  A,  B,  and  C  represent  the  generalized  mass,  added 
mass,  linear  damping,  and  hydrostatic  restoration  coefficients  in  each  equation. 
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Because  of  the  lateral  symmetry  of  the  ship,  several  terms  can  be  set  to  zero. 
The  result  is  that  the  generalized  mass  matrix  has  the  form; 


The  added  mass  matrix  has  the  form: 

Aii  0  Ai3  0  Ai5  0 

0  A22  0  A24  0  A26 

A31  0  A33  0  A35  0 

0  A42  0  A44  0  A46 

A51  0  A53  0  A55  0 

0  A62  0  A64  0  Aee 

The  damping  coefficient  matrix  has  the  same  shape  as  the  added  mass 
matrix.  The  damping  coefficients  are  non-zero  where  the  added  mass 
coefficients  are  non-zero,  and  zero  everywhere  else. 

Almost  all  of  the  restoring  coefficients  for  ships  on  the  ocean  surface  are 
zero.  The  only  non-zero  terms  are  C33,  C44,  C55,  C35,  and  Css-  Expanding  the 
equations  with  all  coefficients  in  them  shows  a  very  interesting  phenomenon. 
Because  of  the  symmetrical  zeroes  in  each  matrix,  the  whole  system  divides  into 
two  sets  of  three  coupled  equations;  one  set  consisting  of  surge,  heave,  and 
pitch,  and  the  other  set  describing  sway,  roll,  and  yaw.  This  can  further  be 
simplified,  as  Salveson,  et.al  point  out[4],  by  noting  that  the  surge  force  and  the 
surge  response,  are  very  small  compared  to  the  other  forces  and  motions.  The 
equations  of  motion  for  each  fixed  frequency  are  then: 
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i^M  +  )^3  +  ^33^3  +  ^33^3  ^35^5  ^35^5  “*"  ^35^5  ^3  (0 
^53^3  “*"  ^53^3  “*"  ^53^3  “*"  (^55  ^55  )^5  ^55^5  ^55^5  ^  ^5  (0 


(3.14) 


and 

(M  +  42)?72  +  ^22%  +  (^4  -Mz^)ij^  +  524/74  +  4g/7g  +  526/7g  =  ^2(0 

(^42  ~  ^^c)V2  "*"  ^42^2  ~*~  (^44  ~*~  ^44)^4  ~*~  ^AaVa  ~*~  ^AaHa  ~*~  (^46  ~  ^  46)^6  ~*~  ^46^6  ^  ^4(0  (3.15) 

(^62)^2  ^62^2  (^64  ~  ^m)Va  “*"  ^64^4  “*"  (^66  ^66)^6  "*"  ^66^6  ~ 

These  five  differential  equations  are  the  basis  of  linearized  seakeeping 
theory  in  the  time  domain.  Because  of  the  assumption  of  linearity,  the  added 
mass  and  damping  coefficients  may  be  considered  functions  of  frequency  and 
forward  speed  only.  The  system  of  equations  represents  a  linear  time 
independent  system,  and  may  be  worked  in  complex  space  by  taking  the  Fourier 
transform  of  each  term.  This,  in  effect,  reduces  the  coupled  differential  equations 
to  algebraic  equations.  The  solutions  to  the  linear  equations  are  in  fact  the 
transfer  functions  for  the  responses,  and  the  response  amplitude  operator  used 
in  spectral  analysis  is  the  magnitude  of  the  transfer  function.  Principles  of  Naval 
Architecture  outlines  the  procedure  for  pitch  and  heave;  the  solutions  provided 
are[5]: 


.  F,S-F,Q 

ri3  =— - ^ 

PS-QR 

.  PP-F.R 
Vi  =— - ^ 

PS-QR 

where 

P  —  C22  —  Q)  (TIT  +  ^4^2 )  +  icoB^^ 
Q  —  C35  —  co^  +  icoB^^ 

R  =  C52  —  CO  A^2  +  icoB^2 
S  —  C55  —  co^  (/jj  +  ^55)  +  icoB^^ 


(3.16) 
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It  should  be  noted  that  the  transfer  functions  are  complex  operators,  with 
only  the  real  part  of  the  product  having  physical  meaning.  An  equivalent 
expression  which  is  used  in  Maxsurf  and  later  in  the  computer  programming  is 
the  modulus  and  phase  of  the  transfer  function //(<») ,  as  shown  below; 

H{Q))  =  ^  =  \H{Q))\e“^^‘"^  (3.17) 

A 

However,  the  transfer  functions  are  only  one  third  of  the  picture.  The 
other  two  parts  are  the  hydrodynamic  coefficients  in  the  equations  and  the 
excitation  forces.  These  values  are  difficult  to  compute,  and  that  difficulty  is 
compounded  by  the  fact  that  many  of  the  coefficients  depend  on  the  frequency  of 
the  incident  waves.  Additionally,  the  excitation  forces  depend  on  both  the 
frequency  and  amplitude  of  the  incident  waves.  The  amplitude  of  the  waves  is 
treated  by  non-dimensionalizing  the  response.  Forces  and  displacements  are 
considered  proportional  to  the  amplitude  of  the  wave,  so  the  non-dimensional 
response  \s  rj!  A  and  the  excitation  for  is  per  unit  amplitude.  Likewise,  the  non- 
dimensional  angular  displacements  arerj/kA . 

As  mentioned,  the  added  mass  and  damping  coefficients  are  dependent 
on  both  the  frequency  and  forward  speed  of  the  ship  in  a  seakeeping  scenario. 
This  thesis  treats  the  added  inertia  and  damping  coefficients  for  roll  as  constants. 
This  simplifies  calculations  and  provides  a  basis  for  comparison  with  the  results 
for  Maxsurf. 

For  pitch  and  heave  motion,  the  Sea  Keeper  module  of  Maxsurf  was  used 
to  obtain  both  the  hydrodynamic  coefficients  and  the  exciting  forces.  The 
resultant  transfer  functions  and  the  discretized  wave  spectrumwere  transformed 
into  time  series  for  both  pitch  and  heave.  The  wave  heights  and  phases  at  many 
discrete  frequencies  were  determined  as  discussed  in  section  3.3.  The  response 
at  time  t  becomes: 


#  freq 

7(0=  X  ^mRAOicoJco^icoJ  +  y/J  (3.18) 

m=\ 
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Examples  of  the  pitch  response  and  heave  response  of  the  ONR 
tumblehome  hull  are  shown  in  the  figures  below.  The  pitch  and  heave  responses 
were  used  as  inputs  to  the  righting  arm  function  discussed  in  section  5.4. 


Figure  3-3:  Heave  Response  in  Sea  State  Eight  Time  Simuiation  from  Sea  Keeper 
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Figure  3-4:  Pitch  Response  in  Sea  State  Eight  Time  Simulation  from  Sea  Keeper 
The  mean  period  for  pitch  is  17.3  seconds,  with  a  root  mean  square 
average  amplitude  of  1 .9  degrees.  This  translates  to  ±5.2  meters  of  pitch  over 
the  length  of  the  ship.  The  heave  period  is  21  seconds,  with  an  average 
amplitude  of  2.4  meters. 

3.4  Treatment  of  Roll  in  MAXSURF 

Maxsurf  uses  the  same  linear  sea  keeping  equations  as  described  in 
Section  3.3  with  three  important  differences. [4]  First,  roll  is  assumed  to  be 
uncoupled  from  sway  and  yaw.  In  other  words,  the  roll  degree  of  freedom  is 
simulated  as  its  own  independent  mass-spring-dashpot  system.  The  reasoning 
for  decoupling  roll  from  sway  and  yaw  is  due  to  the  large  force  and  moment  that 
the  rudder  places  on  a  ship.  Since  this  is  often  ignored  in  seakeeping,  the 
relatively  smaller  wave  force  and  moment  are  ignored.  The  roll  equation  with  this 
assumption  becomes 
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(-^44  ^44)^4  ^AaHa  ^A^ 


(3.19) 


The  second  assumption  that  Maxsurf  makes  deals  with  the  added  mass 
and  damping  coefficients.  In  normal  seakeeping  theory,  these  coefficients  are 
functions  of  excitation  frequency  and  forward  speed.  This  requires  calculations 
in  the  frequency  domain  with  an  inverse  Fourier  transform  (IFFT)  to  obtain  a  time 
domain  solution.  For  Sea  Keeper,  they  are  assumed  to  be  constants.  The 
equation  is  then  an  ordinary  linear  differential  equation  in  the  time  domain  with 
constant  coefficients,  with  the  solution 

n,  =  I —  cos((i)J  +  a)  (3.20) 

V(c„-(/44+4.4)®;f+ii>; 


tan(«)  = 


^44®. 


^"44  (.^44  +  ^44)^^, 


(3.21) 


The  mass  moment  of  inertia  is  calculated  with  a  user-defined  roll  gyradius, 
as  described  in  equation  (2.1)  .  The  added  mass  is  set  to  0.3  times  the  moment 
of  inertia  of  the  ship,  and  the  damping  coefficient  is  set  by  the  user  through  the 
input  of  a  damping  ratio,  P44.  For  this  research,  the  Maxsurf  default  value  of 
0.075  was  used  for  all  calculations.  The  damping  coefficient  is  calculated 
according  to  the  equation 


(3.22) 

Normally,  the  damping  ratio  is  a  function  of  frequency,  making  the 
equation  true  for  each  frequency  separately.  However,  the  frequency 
dependence  is  neglected  by  making  P44  a  constant.  This  permits  a  direct  time 
domain  solution. 
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Finally,  Sea  Keeper  does  not  calculate  the  excitation  moment  using  a 
rigorous  method.  It  assumes  that  the  moment  may  be  approximated  as  the 
product  of  the  wave  slope  and  the  hydrostatic  righting  moment,  in  phase  with  the 
wave  slope.  Therefore: 


F4  =kAC^^ 

(3.23) 

where 

c^=  pgyGM, 

(3.24) 

Using  all  of  the  values  from  this  section,  the  roll  response  amplitude 
operator  for  roll  becomes 


Il±  =  RAO  = 
kA 


C 


44 


^/(C 


44  “  (^4  +  ^44)®^)^  + 


44 


(3.25) 


The  roll  response  amplitude  operator  may  be  used  in  conjunction  with  the 
wave  height  spectrum  to  produce  a  time  series  simulation  just  as  in  pitch  and 
heave.  An  example  of  the  time  simulation  for  roll  is  shown  below: 
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Figure  3-5:  Roll  Response  Time  Simulation  from  Sea  Keeper 


Although  Sea  Keeper’s  treatment  of  roll  does  not  give  the  most  accurate 
results,  it  does  give  a  basis  for  comparison  of  the  results  of  this  research. 
Therefore,  the  new  programs  for  making  time  simulations  of  roll  angle  utilized  the 
same  assumptions  for  moment  of  inertia,  added  inertia,  and  damping  that  Sea 
Keeper  uses.  The  modifications  used  in  this  research  are  to  the  restoring 
coefficient  C44  and  the  excitation  moments.  The  non-linearities  of  the  restoring 
coefficients  have  been  discussed  in  section  2.3,  and  the  treatment  of  the 
excitation  forces  will  be  discussed  in  the  next  section. 

4  Roll  Excitation  Moment 

Many  programs  provide  a  spectrum  for  the  excitation  roll  moment  due  to 
incident  waves.  Sea  Keeper’s  gives  an  order  of  magnitude  estimation  of  the 
excitation  roll  moment.  A  better  estimate  can  be  obtained  using  the  linearized 
potential  due  to  incident  waves.  The  water  particles  exert  a  dynamic  pressure  on 
the  hull  surface,  which  can  be  integrated  over  the  surface  of  the  ship  to  give  a 
good  approximation  of  the  excitation  moment.  Typically,  the  moment  is 
calculated  in  two  parts.  First,  the  moment  due  to  the  incident  potential  alone, 
called  the  Froude-Krylov  moment,  is  calculated.  The  second  part  calculated  is 
the  moment  due  to  diffraction  of  the  wave  as  it  interacts  with  the  ship.  Both  the 
Froude-Krylov  moment  and  the  Diffraction  moment  will  be  discussed  in  this 
section. 

4.1  Froude-Krylov  Moment 

Linear  plane  progressive  waves  are  assumed  to  impact  the  ship  at  an 
angle  p  to  the  stern  of  the  ship.  As  discussed  previously,  the  two  dimensional 
incident  wave  potential  is  given  by  equation(3.7) .  The  moment  resulting  from  the 
incident  potential  is  given  by  the  Froude-Krylov  hypothesis,  which  is  expressed 
as: 
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In  complex  space,  the  Froude-Krylov  hypothesis  is  written: 

O.n.dS  (4.2) 

C’ 

In  this  case,  the  surface  C  is  the  wetted  surface  of  the  ship’s  hull  in  still 
water,  and  the  normal  vector  component,  r\4,  is  the  cross  product  of  the  hull’s 
two-dimensional  normal  vector  with  the  position  vector  of  the  point  in  question, 
shown  in  Figure  4-1: 


n4  =  I  r  X  n  | 

Figure  4-1 :  Normal  Vector  for  Moments 

The  method  of  numerical  analysis  becomes  quite  clear.  For  a  given 
section,  the  normal  and  incident  potential  are  known  at  each  point  in  the  table  of 
offsets.  The  complex  sectional  Froude-Krylov  force  is  expressed  as: 


ft  points 


#  po  int  s 


(4.3) 


And  the  overall  moment  can  be  calculated  using  a  trapezoidal  rule  along  the 
length  of  the  ship: 
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(4.4) 
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4.2  Diffraction  Moment 


The  excitation  moment  due  to  the  diffraction  of  waves  around  the  hull  of 
the  ship  is  more  complicated.  This  thesis  uses  the  calculation  of  the  diffraction 
moment  used  by  Milgram  [6].  Like  the  Froude-Krylov  moment,  the  diffraction 
moment  may  calculated  as 


<=-4 


ao 


dt 


^n^dS 


(4.5) 


The  problem  is  that  the  diffraction  potential  is  not  known,  and  more  than 
likely  can  not  be  found  as  an  explicit  function  of  position  and  time.  Therefore, 
numerical  methods  must  be  used  to  determine  the  diffraction  potential  before 
numerical  integration  can  be  implemented. 

4.2. 1  Diffraction  Boundary  Value  Problem 


Just  as  in  the  case  of  plane  progressive  waves,  the  diffraction  potential 
must  satisfy  a  fluid  flow  boundary  value  problem.  The  boundary  value  problem 
can  be  stated  as  follows:  the  total  potential  must  satisfy  mass  conservation 
inside  the  fluid,  with  no  flux  along  the  hull  of  the  ship  or  bottom  boundary  of  the 
domain.  The  potential  of  the  sides  and  free  surface  reflect  the  fact  that  the 
diffraction  will  produce  waves  radiating  outward.  Therefore,  the  diffraction 
potential  satisfies  the  following  equations: 
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4.2.2  Numerical  Determination  of  Diffraction  Potential 

Once  the  boundary  conditions  are  known,  the  potential  may  be 
determined  using  a  constant  panel  method  using  Green’s  theorem.  The  methd 
for  determining  the  diffraction  potential  is  modeled  on  that  used  by  Milgram[6].  It 
is  known  that  the  potential  obeys  Laplace’s  equation  in  the  field.  Another 
function  that  satisfies  LaPlace  is  required  for  Green’s  theorem.  In  this  case,  the 
Rankine  Green’s  function  was  used.  The  Rankine  Green’s  function  takes  the 
form 

G{y,  z,r/,g)  =  -ln(yl(y-r/f  +(z-gf)  (4.11) 


This  function  therefore  depends  on  the  distance  from  the  field  point  (y,z) 
and  the  source  point  (ri,^).  Because  both  the  diffraction  potential  and  Green’s 
function  satisfy  LaPlace’s  equation.  Green’s  theorem  states  that  the  following 
equation  is  true  along  the  boundary  of  the  domain: 


dG 

dn(rj,g) 
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d^{ri,g) 
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)ds  =  -7irt>{y,z) 


(4.12) 


To  solve  this  numerically,  the  domain  must  be  divided  into  several  panel 
line  elements.  Along  each  panel,  the  values  of  O  and  6ct>/6n  in  the  equation 
above  may  be  assumed  to  be  constant.  This  means  that  the  panels  must  be 
small  enough  that  the  potential  does  not  change  significantly  across  the  panel. 
This  becomes  important  is  selecting  the  domain  size,  since  more  panels 
significantly  increases  the  computation  time  for  each  section.  The  result 
becomes  a  summation  throughout  the  domain,  namely; 
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The  integrals  in  the  above  equation  can  be  computed  numerically  as  well. 
The  boundary  conditions  provide  relationships  to  eliminate  the  derivative  of  the 
diffraction  potential  at  each  point.  This  reduces  the  system  to  a  set  of  N  linear 
equations  for  the  diffraction  potential. 

For  the  problem,  a  two  dimensional  domain  chosen  for  this  problem  is  a 
box  of  ocean  fluid  with  the  hull  in  the  center  of  the  top,  as  shown  in  Figure  4-2: 


Figure  4-2:  Domain  for  Diffraction  Boundary  Vaiue  Probiem 

The  half-width  of  the  domain  W  was  set  to  ten  ship  beams.  This  width 
was  determined  to  be  acceptable  for  convergence  by  Milgram  [6].  The  depth  FI 
was  set  to  one  half  the  wavelength  of  the  incident  wave  or  three  times  the  draft  T, 
whichever  was  greater.  This  brings  the  domain  to  where  the  incident  potential  is 
near  zero  for  long  waves,  and  covers  the  bottom  of  the  hull  for  short  waves.  The 
number  of  panels  along  the  free  surface  balanced  computation  time  against 
panel  length.  At  55  panels  on  either  side  of  the  ship,  the  panel  length  is 
approximately  1.5  meters.  With  a  cutoff  frequency  of  2  radians  per  second,  this 
means  that  each  panel  is  less  than  one  tenth  of  the  wavelength  even  at  the 
shortest  wavelength. 

5  Matlab  Code  Description 

The  objective  for  the  computer  programs  is  to  obtain  a  reasonable 
estimate  to  the  non-linear  roll  response  using  a  desktop  computer  without 
requiring  unreasonable  computation  time.  For  the  most  part,  the  objective  is 
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accomplished  well.  The  programming  language  used  is  Matlab,  release  14.  The 
computer  used  was  a  Dell  Optiplex  GX620,  with  a  3.2  GHz  Pentium  4  processor 
and  1GB  of  memory. 

5.1  Files  Required  from  Maxsurf 

The  computer  codes  written  for  this  research  require  three  files  containing 
output  from  Maxsurf.  First,  a  table  of  offsets  for  the  hull  in  question  must  be 
generated  for  use  in  the  force  calculator  program.  This  file  is  in  Microsoft  Excel 
format,  with  four  columns.  The  first  column  is  the  station  number  of  the  offset 
point.  The  program  is  designed  to  handle  an  arbitrary  number  of  stations.  The 
next  three  columns  are  x-coordinates,  z-coordinates,  and  y-coordinates, 
respectively.  The  reason  for  the  switch  in  y-  and  z-coordinates  is  that  the  offsets 
may  be  cut  and  pasted  directly  from  a  table  in  Maxsurf.  The  title  of  the  file  must 
be  ‘offsets.xls.’  Generation  of  this  file  in  Maxsurf  took  approximately  four  hours. 
Most  of  this  time  was  generating  offsets  from  the  imported  surface. 

The  next  file  required  for  the  programs  is  used  in  the  linear  response 
generator.  This  file  contains  the  response  amplitude  operators  for  pitch,  heave, 
and  roll.  Each  RAO  is  defined  as  a  500  element  row  vector  for  amplitude  and 
phase.  Additionally,  the  corresponding  frequencies,  encounter  frequencies,  and 
wave  amplitudes  calculated  as  described  in  section  3.2  are  included  in  the  file. 
The  file  therefore  contains  nine  row  vectors  with  the  following  names:  HRAO 
(heave  response);  Hph  (heave  phase);  RRAO  (roll  response);  Rph  (roll  phase); 
PRAO  (pitch  response);  Pph  (pitch  phase);  freqs  (wave  frequencies);  efreq 
(encounter  frequencies);  and  A  (wave  amplitudes).  The  name  of  the  file  is 
‘raofunctions.mat.’  This  file  required  about  one  hour  to  generate  once  the  hull 
was  saved  as  a  Maxsurf  file. 

A  complication  to  generating  this  file  was  the  length  of  the  vectors.  The 
simulation  included  frequencies  up  to  2  radians  per  second,  and  500  discrete 
frequencies  were  desired  to  provide  adequate  randomness  for  a  24  hour  period. 
Since  Maxsurf  could  only  produce  132  frequencies  in  this  range,  Matlab’s  interpi 
function  was  used  to  interpolate  values  for  the  500  elements  of  each  variable. 
This  did  not  significantly  increase  the  preparation  time  for  the  file. 
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Finally,  the  righting  moment  lookup  table  is  required.  This  is  a  three- 
dimensional  lookup  table  with  reference  values  for  each  index.  The  method  for 
generating  the  lookup  table  utilized  Maxsurf’s  Hydromax  module.  At  each  draft, 
pitch  was  fixed  at  values  from  20m  by  the  stern  to  20m  by  the  bow  (+/-  .133 
radians.  This  ensures  all  reasonable  values  of  pitch  angle  in  the  time  simulation 
can  be  accommodated.  The  data  was  collated  in  a  spreadsheet  table  in  the 
format  shown  below: 


Displ(kg)  LCG(m)  VCG(m)  T(m) 

1.6E+07  -10.8  7.5  9 


Roll(degrees) 


Trim(m) 

0 

10 

20 

30 

40 

20 

0 

3143262 

6111013 

8727739 

10754106 

15 

0 

3079440 

6126968 

8951118 

11216820 

10 

0 

2584815 

5552565 

8647961 

11312554 

5 

0 

1850855 

4244202 

7690622 

11137042 

0 

0 

1595565 

3685754 

6797106 

10770062 

-5 

0 

1771077 

4004867 

7068352 

10227570 

-10 

0 

2185924 

4659049 

7195997 

9174497 

-15 

0 

2313569 

4531404 

6398214 

7770400 

-20 

0 

2058278 

3893178 

5408964 

6462037 

Displ(kg)  LCG(m)  VCG(m)  T(m) 
1.5E+07  -10.625  7.5 

8.75 

Roll(degrees) 

Trim(m) 

0 

10 

20 

30 

AH 

20 

0 

3208086 

6246513 

8914777 

10966101 

15 

0 

3146392 

6261937 

9146129 

11459652 

10 

0 

2544876 

5537033 

8745118 

11521346 

5 

0 

1804548 

4148919 

7557510 

11197453 

0 

0 

1542349 

3578249 

6601253 

10611360 

-5 

0 

1712007 

3855872 

6755488 

9948150 

-10 

0 

2066748 

4395694 

6940570 

8961047 

-15 

0 

2205559 

4364847 

6246513 

7650051 

-20 

0 

2005054 

3825025 

5336527 

6385324 

Table  5-1:  Righting  Moment  vs.  Roll,  Pitch,  &  Draft 


The  table  contains  restoring  moment  divided  by  the  acceleration  due  to 
gravity  (mass  times  righting  arm)  for  each  value  of  heave,  pitch,  and  roll.  The 
two-dimensional  table  of  righting  moment  as  a  function  of  roll  and  pitch  are  a  two- 
dimensional  table  for  the  i‘^  set  of  the  variable  RM.  The  index  vectors  contain  the 
sequential  values  of  trim,  roll,  and  draft.  The  titles  of  these  vectors  are  drafts, 
pitchs,  and  rolls.  The  pitch  is  in  radians  from  high  to  low,  while  the  roll  is  in 
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radians  from  0  to  1 .74  radians.  Drafts  are  in  0.25  meter  increments  from  two 
meters  to  nine  meters.  This  file  is  used  in  the  non-linear  response  generator. 
The  entire  file  for  nine  pitches,  eleven  roll  angles,  and  29  drafts  took 
approximately  eight  hours  to  compile. 

5.2  Excitation  Force  Calculator 

The  Matlab  program  ForceFinder  calculates  the  linear  excitation  moment 
of  a  ship  from  offsets  provided  in  an  external  file  using  the  theory  presented  in 
Chapter  Four.  The  program  is  designed  to  calculate  the  Froude-Krylov  and 
Diffraction  excitation  moments  for  a  user  defined  number  of  different  frequencies. 
As  mentioned,  the  program  requires  a  table  of  offsets  with  points  designated  in 
meters  from  the  origin.  At  each  frequency,  the  program  determines  sectional 
moments  per  unit  length  of  the  ship  at  each  section.  These  sectional  moments 
are  then  integrated  over  the  length  of  the  ship  to  obtain  the  three  dimensional 
excitation  force.  A  block  diagram  of  the  program  is  shown  in  figure  5.1 . 


Figure  5-1:  ForceFinder  Biock  Diagram 


ForceFinder  uses  several  subroutines  to  determine  the  sectional 
moments.  The  main  subroutine  is  ‘findfk’.  This  subroutine  takes  a  section’s 
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offsets,  determines  the  sectional  excitation  moments  for  the  frequency  passed  to 
it.  First,  it  generates  the  domain  shown  in  figure  4-2  using  the  subroutine 
‘setpans’.  Next,  the  subroutine  ‘inficoef  calculates  the  integrals  of  Green’s 
function  and  its  derivative  at  each  point  due  to  every  other  point.  The  subroutine 
‘findmatrix’  uses  these  values  and  the  boundary  conditions  to  calculate  the 
equations  for  determination  of  the  diffraction  potential.  With  the  diffraction 
potential  known,  the  subroutine  calculates  both  the  Froude-Krylov  and  diffraction 
moments  through  numerical  integration. 

Although  this  process  is  straightforward,  a  block  diagram  of  findfk  is 
presented  in  figure  5-2.  The  code  for  ForceFinder,  findfk,  and  findmatrix  are 
included  in  section  1  of  Appendix  A. 


Figure  5-2:  findfk  Block  Diagram 

For  the  research  program,  the  program  was  set  to  generate  excitation 
moments  at  500  frequencies  equally  divided  from  zero  to  2  radians  per  second. 
This  matches  the  number  of  wave  amplitudes  generated  for  the  linear  response 
spectrum.  In  this  way,  the  random  phasing  in  the  linear  response  may  be  applied 
to  the  excitation  moments  as  well.  At  each  frequency,  the  moment  and 
amplitude  for  a  one  meter  wave  was  calculated.  The  result  is  a  complex  number, 
which  may  be  interpreted  as  a  modulus  of  moment  and  a  phase  relative  to  the 
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wave  height  amidships.  The  output  of  the  program  was  saved  as  a  Matlab  data 
file  titled  ‘forcedataSOO.mat’  for  use  in  the  non-linear  time  simulations.  The 
excitation  moment  per  unit  wave  amplitude  for  Sea  State  Eight  is  shown  in  figure 
5-3: 
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Figure  5-3:  Excitation  Moments  for  Seakeeping  Problem 


ForceFinder  has  the  largest  run  time  of  the  three  programs  described. 
Because  of  the  large  number  of  panels  in  the  diffraction  boundary  value  problem, 
the  program  requires  approximately  50  seconds  to  find  the  total  excitation  force 
for  each  frequency.  To  match  the  500  frequencies  used  in  the  spectrum 
response  calculator,  the  program  run  time  was  6.9  hours.  The  author 
acknowledges  there  are  many  programs  available  to  evaluate  the  excitation  force 
that  are  faster.  These  may  be  used  as  long  as  the  output  is  in  the  form  of  a 
Matlab  vector  labeled  F4  containing  complex  force  magnitudes  and  is  saved  as 
‘forcedata500’. 
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5.3  Linear  Response  Generator 


The  name  of  the  linear  response  program  is  LinearResponse.  Although  it 
does  not  need  any  subroutines  to  operate,  it  does  need  the  file  ‘raofu notions. mat’ 
to  generate  a  time  simulation  of  heave,  pitch,  and  roll.  Although  only  the  heave 
and  pitch  vs.  time  are  required  for  use  in  the  non-linear  program,  the  roll  was 
calculated  for  comparison  to  the  time  simulation  generated  in  the  non-linear  roll 
simulation. 

The  linear  response  program  and  force  program  should  be  run  with  the 
same  set  of  frequencies  in  their  spectra.  The  linear  response  program  calculates 
the  wave  amplitudes  and  random  phases  for  each  frequency  in  the  simulation. 
The  moments  found  in  ForceFinder  have  a  phase  relative  to  the  phase  of  their 
respective  discrete  wave  system.  If  the  frequencies  are  not  matched  prior  to 
running  the  two  programs,  a  certain  amount  of  interpolation  must  happen  to 
synchronize  the  frequencies  before  the  final  program  can  be  run.  This  extra  work 
is  easily  avoidable  by  using  the  same  set  of  frequencies  for  the  separate 
programs. 

The  linear  response  program  is  quite  simple.  For  each  time,  the  response 
may  be  expressed  as  a  sum  of  the  linear  response  at  each  frequency.  This  may 
be  written  as 

i^freq 

=  S  4  \H{o}i)\ cos{o)f  +  y/^  +  «,.)  j=3,4,5  (5.1) 

i=\ 

The  file  raofunctions.mat  described  in  section  5.1  contains  the  wave 
amplitudes,  transfer  function  moduli,  and  transfer  function  phase  angles  for  these 
calculations.  Therefore,  the  program  generates  a  random  phase  angle  for  each 
frequency,  and  performs  the  summation  of  discrete  responses  at  each  time  to 
produce  a  wave  height,  pitch  angle,  roll  angle,  and  heave. 

The  output  of  this  program  is  the  time  series  simulations  of  wave  height, 
pitch,  heave  and  roll.  Also  in  the  data  set  will  be  the  frequencies,  encounter 
frequencies,  and  wave  amplitudes  for  the  sea  simulation.  The  output  was  saved 
as  a  file  titled  ‘linresplOO.mat’  for  use  in  the  non-linear  response  simulation 
program.  Sample  graphs  of  the  time  simulations  were  shown  previously  in 
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Figure  3-3,  Figure  3-4,  and  Figure  3-5.  Run  time  for  this  program  is 
approximately  40  seconds.  The  code  is  included  in  Appendix  A. 

5.4  Non-Linear  Response  Generator 

The  name  of  the  non-linear  response  program  is  ‘rollintegration.m’.  It 
uses  one  subroutine,  a  three  dimensional  lookup  and  interpolation  routine  to 
determine  the  righting  moment  as  a  function  of  pitch,  roll,  and  heave.  For  this 
program  to  work  properly,  it  requires  the  output  of  ForceFinder.m  saved  as  a  file 
titled  ‘forcedata500.mat’,  the  output  of  LinearResponse.m  saved  as 
‘linresplOO.mat’,  and  the  righting  moment  lookup  table  saved  as 
‘Rmomentdata.mat’. 

The  program  uses  a  forward  Euler  scheme  to  integrate  the  initial  value 
problem  defined  by  equation  (3.19).  Flere,  however,  C44  at  each  time  step  is  the 
non-linear  function  C44 [773(0, 774 (r), //s (r)] .  Since  this  equation  is  a  second  order 
differential  equation,  it  must  be  treated  as  two  first  order  equations  to  evaluate  it 
numerically. [7]  Therefore,  using  dummy  variables  qi  and  q2,  equation  (3.19) 
becomes 

(5.2) 


At  each  time  step,  the  excitation  moment  is  determined  by  the  equation 


H-freq 

M{t)=  '^\F^+H^\cos{Q}J  +  y/ (5.3) 

n=\ 


Initial  conditions  for  both  roll  angle  and  its  first  derivative  are  chosen  to 
match  the  output  from  ‘linresplOO.mat’.  The  second  time  derivative  is  evaluated 
for  the  present  time  step  and  then  integrated  to  give  the  first  derivative  at  the  next 
time  step.  The  first  time  derivative  at  the  present  time  step  is  integrated  to 
calculate  the  roll  angle  at  the  next  time  step.  At  any  time  step,  if  the  non-linear 
roll  angle  grows  to  greater  than  ninety  degrees,  the  counter  for  catastrophic  roll 
events  is  advanced  by  one  and  the  roll  angle  and  speed  are  reset  to  the  values 
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of  the  linear  integration  at  for  that  time.  The  block  diagram  for  this  program  is 
shown  in  Figure  5-4. 


Figure  5-4:  Block  Diagram  for  Non-Linear  Roll  Integration 


An  additional  capability  written  into  the  program  is  to  conduct  the 
integration  using  the  same  linear  restoring  coefficient  that  Maxsurfs  Sea  Keeper 
module  uses,  using  equation(3.24).  The  linear  and  non-linear  integration  occur 
in  the  same  loop,  so  little  extra  computation  time  is  required  for  this  addition.  The 
programs  as  written  in  Appendix  A  are  designed  to  use  500  frequencies  up  to  a 
value  of  2  radians  per  second  to  integrate  the  roll  equation  for  one  hour  using  a 
time  step  of  0.01  seconds.  The  run  time  for  rollintegration  is  approximately  one 
minute.  The  small  time  step  is  required  for  numerical  convergence.  For  such  a 
small  time  step,  the  Forward  Euler  method  gives  very  nearly  the  same  results  as 
more  complicated  integration  schemes  such  as  Fourth  Order  Runge-Kutta. 
Convergence  with  a  larger  time  step  for  implicit  integration  rules  remains  to  be 
investigated  for  solving  the  roll  differential  equation. 
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6  Simulation  Results 

6.1  Sea  Keeper  and  Linear  Integration  Simulations 


The  first  comparison  that  can  be  made  using  the  output  of  the  integration 
program  is  between  the  Sea  Keeper  module  of  Maxsurf  and  the  integration  using 
the  linear  forces  and  a  linear  restoring  coefficient.  As  noted  in  section  3.4,  Sea 
Keeper  uses  an  approximation  of  wave  slope  times  the  hydrostatic  restoring 
coefficient  for  its  excitation  force.  To  compare  the  two  excitation  forces,  the 
output  from  the  ForceFinder  program  was  normalized  by  the  value  of  the  Sea 
Keeper  excitation  force  (kAC44).  This  comparison  is  shown  in  Figure  6-1 .  The 
figure  shows  the  normalized  modulus  of  the  excitation  force  vs  frequency,  along 
with  the  constant  value  of  kAC44  that  Sea  Keeper  uses  as  its  excitation  force.  The 
largest  discrete  wave  heights  will  be  near  the  modal  frequency,  which  is  also 
shown  in  the  figure. 
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Figure  6-1:  Modulus  of  Non-dimensional  Excitation  Force  vs.  Frequency 


The  phase  angles  of  the  excitation  forces  are  not  shown.  Since  the 
calculated  excitation  moment  is  computed  as  a  complex  number,  the  phase 
varies  with  the  frequency  of  the  wave.  On  the  other  hand,  the  Sea  Keeper 
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excitation  moments  do  not  include  a  phase  shift,  assuming  that  the  excitation 
force  always  acts  in  phase  with  the  wave  height.  While  this  assumption  may  be 
valid  for  longer  wavelengths,  it  does  not  hold  true  for  shorter  waves.  This 
difference  does  not  significantly  affect  the  amplitude  of  the  response. 

The  modulus  of  the  calculated  force  is  larger  in  the  area  of  the  largest 
discrete  wave  heights.  Therefore,  the  overall  calculated  force  should  be  larger 
than  the  estimations  used  by  Sea  Keeper.  The  linear  integration  uses  the  same 
value  of  C44  as  the  Sea  Keeper  program.  Therefore,  the  only  difference  of  the 
time  series  for  roll  between  the  integration  scheme  and  Sea  Keeper  is  the 
excitation  forces.  The  larger  forces  in  the  integration  should  result  in  a  larger 
average  amplitude  for  roll  response.  This  is  in  fact  the  case.  The  average 
amplitude  was  determined  by  calculating  the  standard  deviation  of  the  time 
series  of  the  two  different  methods.  The  amplitude  of  the  Sea  Keeper  time  series 
was  8.8  degrees,  while  the  linear  integration  average  amplitude  was  11  degrees. 
A  plot  of  the  first  five  minutes  is  shown  in  Figure  6-2  for  comparison  of  the  two. 
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Figure  6-2:  Roll  Angle  vs  Time  from  Linear  Integration  and  Response  Spectrum 


The  period  can  be  approximated  by  dividing  the  time  period  by  the  number 
of  zero  crossings  in  the  simulation.  If  this  is  done,  the  roll  period  for  the  linear 
integration  is  17.3  seconds,  while  the  period  of  the  Sea  Keeper  simulation  is  14.1 
seconds.  Both  of  these  agree  with  the  natural  period  estimated  section  2.2.  The 
difference  in  them  reflects  the  concentration  of  roll  excitation  at  different 
frequencies  due  to  the  different  functions  of  the  roll  excitation  moment.  The 
results  of  the  two  different  methods  are  rather  close,  showing  that  the  force 
calculation  program  appears  to  give  results  that  are  valid.  The  roll  period  is  long 
because  of  the  forward  speed  effects  in  stern  quartering  seas.  The  forward 
speed  concentrates  the  sea  state’s  excitation  forces  at  lower  frequencies, 
causing  a  longer  roll  period. 

6.2  Linear  Integration  and  Non-Linear  Integration 

The  linear  time  simulation  and  the  non-linear  time  simulation  have  only 
one  difference.  The  restoring  force  in  the  non-linear  case  is  a  function  of  pitch 
and  heave,  as  well  as  a  function  of  roll.  As  noted  in  section  2.3,  the  restoring 
coefficient  could  change  by  up  to  a  factor  of  three  due  to  pitch  about  its  still 
waterline.  The  response  should  have  approximately  the  same  frequency 
behavior,  but  there  should  be  periods  where  the  amplitude  of  the  linear  model  is 
larger  than  that  of  the  non-linear  model,  and  periods  where  the  reverse  is  true. 

Surprisingly,  very  few  cases  could  be  found  where  the  linear  response  had 
larger  amplitude.  In  Figure  6-3  below,  three  such  occurrences  happen  in  five 
minutes.  As  you  can  see,  the  larger  linear  response  is  not  much  larger  than  the 
non-linear  response.  Conversely,  Figure  6-4  shows  the  much  more  frequent 
case  of  the  larger  non-linear  response.  The  non-linear  response  peaks  tend  to 
be  much  larger.  This  indicates  that  the  non-linear  effect  is  a  net  loss  in  stability 
for  the  ship  as  it  pitches  and  heaves  in  stern  quartering  seas. 
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Figure  6-3:  Linear  Simulation  vs.  Non-linear  Simulation  with  Larger  Rolls  in  the  Linear 

Case 

The  period  of  the  two  integration  schemes  were  very  similar.  While  the 
linear  integration  period  was  17.3  seconds,  the  non-linear  time  series  showed  a 
period  of  17.0  seconds.  The  amplitudes  of  the  responses  were  quite  different, 
however.  As  noted  before,  the  standard  deviation  of  the  linear  series  was  1 1 
degrees.  The  non-linear  roll  had  a  standard  deviation  of  16.9  degrees,  which  is 
an  increase  of  more  than  fifty  percent.  With  larger  roll  amplitude  and  the  same 
period,  this  implies  much  greater  roll  speed  and  acceleration.  Therefore  not  only 
is  the  hull  less  stable  when  treated  in  a  non-linear  fashion,  it  would  be  less 
comfortable  for  crew  as  well. 
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Figure  6-4:Linear  Simulation  vs.  Non-linear  Simulation  with  Larger  Rolls  in  the  Non-linear 

Case 

Examination  of  Figure  6-4  shows  a  phenomenon  that  is  unique  to  the  non¬ 
linear  case.  The  roll  angle  exceeds  ninety  degrees.  This  large  of  a  roll  angle 
exceeds  the  ability  of  the  program  to  determine  the  righting  moment,  and  can  be 
called  a  catastrophic  roll  event.  This  hull  form  may  have  a  positive  righting 
moment  at  roll  angles  of  more  than  ninety  degrees,  but  a  roll  of  more  than  ninety 
degrees  could  be  catastrophic  nonetheless  for  two  reasons.  First,  anything 
inside  the  ship  that  is  not  tied  down  will  undoubtedly  move  towards  the 
furthermost  hull.  Among  other  problems  this  represents,  the  change  in  the  center 
of  gravity  will  cause  an  upsetting  moment  within  the  ship.  The  second 
undesirable  effect  would  be  downflooding.  The  downflood  angles  of  this  hull 
design  are  not  designated,  but  immersing  half  of  the  main  deck  and 
superstructure  should  almost  certainly  cause  down  flooding  even  for  roll  angles 
less  than  90  degrees. 
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Overall,  the  ship  experienced  two  catastrophic  roll  events  in  the  one  hour 
simulation  during  the  non-linear  case.  Neither  linear  simulation  contained  a 
similar  event.  A  comparison  of  average  period  and  amplitude  for  the  three  cases 
is  shown  in  Table  6-1 . 


Simulation 

Period(sec) 

Amplitude(degrees) 

Catastrophic  Events 

Sea  Keeper 

14.1 

8.8 

0 

Linear 

17.3 

11 

0 

Integration 

Non-Linear 

17.1 

16.9 

2 

Integration 

Table  6-1:Comparison  of  Different  Roll  Response  Calculations 


One  interesting  use  for  the  integration  program  was  to  evaluate  the 
influence  of  initial  trim  on  the  stability  of  the  program.  A  0.5  meter  trim  by  the 
stern  was  added  to  the  pitch  response  before  the  time  series  for  roll  was 
generated.  Even  this  small  amount  of  pitch  was  beneficial.  Trimming  the  ship 
reduced  the  mean  amplitude  of  the  roll  by  three  degrees  in  the  non-linear  case 
and  reduced  the  number  of  catastrophic  events  to  one.  Obviously,  in  a  real 
situation,  additional  actions  can  be  taken  to  prevent  such  large  roll  events.  As 
stated  in  section  2.4,  the  high  seas,  speed  of  the  ship,  and  angle  of  incidence  of 
the  incoming  waves  are  all  chosen  to  evoke  a  maximum  roll  response. 
Therefore,  slowing  the  ship  and  maneuvering  for  bow  seas  would  further  improve 
the  ship’s  seakeeping  ability. 

7  Future  Work  and  Conclusion 

7.1  Added  Mass  and  Damping 

As  pointed  out  in  section  3.3,  the  added  mass  and  damping  coefficients 
were  treated  as  constants  instead  of  functions  of  frequency.  An  addition  to  the 
excitation  moment  calculator  program  could  be  the  solution  to  the  radiation 
potential  at  each  frequency,  which  will  give  the  added  mass  and  damping 
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coefficients  for  use  in  the  non-linear  integration  problem.  The  frequency 
dependent  added  mass  and  damping  can  be  calculated  by  a  numerical  problem 
very  similar  to  that  used  to  find  the  diffraction  force[6].  Only  the  boundary 
conditions  on  the  ship  surface  are  different. 

7.2  Directional  Seas 

The  waves  that  are  simulated  in  this  thesis  are  long-crested  waves 
assumed  to  be  from  a  distant  storm.  Simulation  of  multidirectional  random  seas 
may  be  included  in  this  analysis  through  methods  described  in  previous  research 
from  Mr.  Sam  Geiger[8].  The  unidirectional  spectrum  could  be  replaced  with  a 
sum  of  several  unidirectional  seas  with  different  headings,  all  adding  up  to  the 
same  energy  as  the  single  unidirectional  sea  state.  This  inclusion  would  provide 
a  better  estimate  of  the  excitation  forces  involved  in  real  seas. 

7.3  Non-linear  Excitation  Moments 

The  Froude-Krylof  and  Diffraction  excitation  moments  were  calculated 
using  linear  strip  theory  along  the  still  waterline  of  the  ship’s  hull.  In  reality,  the 
moments  will  depend  on  the  actual  waterline  of  the  ship  as  the  waves  travel  in 
space.  The  moment  will  also  depend  on  the  amount  of  the  hull  in  the  water  as 
the  ship  is  displaced  in  roll,  heave,  and  pitch.  Calculation  of  these  effects  is  left 
as  a  topic  of  further  research. 

7.4  Wave  Profile  Effects  on  Righting  Moment 

The  righting  moment  for  the  tumblehome  hull  was  calculated  using  the 
draft  of  the  ship  at  amidships,  assuming  no  effects  for  wave  height  or  wave 
slope.  The  draft  at  other  positions  along  the  length  of  the  ship  is  a  function  of  the 
longitudinal  wave  height  and  the  pitch  of  the  ship  relative  to  the  wave  slope. 
Therefore,  the  righting  moment  may  be  refined  by  taking  these  two  factors  into 
consideration. 

7.5  Conclusion 
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The  non-linear  roll  integration  program  developed  in  this  thesis  can 
provide  an  estimate  of  the  roll  response  for  any  hull,  given  its  offsets  and  the 
righting  moment  data.  Although  the  force  calculation  program  takes  some 
computation  time,  any  available  program  that  determines  the  excitation  force 
may  be  used  to  give  the  roll  integration  program  its  input.  The  program  could  be 
used  in  any  number  of  design  applications,  such  as  determination  of  a  safe 
operating  envelope  for  a  ship,  or  the  effects  of  different  trim  configurations  on 
stability.  Most  importantly,  it  provides  a  non-linear  simulation  in  a  relatively  short 
time  using  only  a  desktop  computer.  This  fact  makes  the  computer  program  a 
valuable  design  tool. 
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Appendix  A:  Matlab  Code  for  Finding  Excitation  Forces 


1.  Main  Program:  ‘ForceFinder.m’  -  This  program  needs  the  following 
subroutines  to  function  :  findfk,  findmatrix,  setpans,  inficoef,  rank2d,  localize. 
Also  needs  a  Microsoft  Excel  file  with  offsets  provided  as  described  in  the 
program.  The  m-files  inficoef,  rank2d,  and  localize  were  provided  by  Professor 
Jerome  Milgram. 

%  Finds  Froude-Krylov  and  Diffraction  Excitation  Moments 
%  for  frequencies  up  to  2.5  rad/s,  unit  amplitude  wave 
%  Requires  the  following  files  to  operate; 

%  offsets.xls;  table  of  offsets  in  MS  Excel  with  the  following  columns: 

%  Col  1 :  station  number  (can  handle  any  number  of  stations) 

%  Col  2:  x-values  (m) 

%  Col  3;  z-values  (m) 

%  Col  4:  y-values  (m) 

%  Subroutines:  findfk,  findmatrix,  setpans,  inficoef,  rank2d,  localize 
%  #  of  frequencies  is  set  by  nfreq 
%  output  is  in  the  variables  F4,  H4 

%  for  the  research,  set  nfreq  to  500,  save  data  as  forcedata500 

clear  all 
close  all 
tic 

timestampb=clock; 

starttime=timestampb(4)*100+timestampb(5); 
fprintf('start  time  =%1 .0f\n', starttime) 

%  initialization  of  variables: 

%  rho:  density  of  seawater  (kg/m''3) 

%  nsta;  #  of  stations  (from  offsets  file) 

%  B;  beam  of  ship  (m) 

%  gr:  acceleration  due  to  gravity 
%  npts:  number  of  offsets 
%  U;  foHA/ard  speed  of  the  ship 
%  nfreq;  number  of  frequencies  in  interval  0:2.5 
%  F4:  Foude  Krylov  excitation  moment 
%  H4:  Diffraction  excitation  moment 

%  wo,  we,  k:  stationary  frequency,  encounter  frequency,  wave  number 
%  T:  draft 

%  beta:  angle  of  incidence  of  waves 
rho=1025; 

A=  xlsread('offsets'); 
nsta=max(A(:,1)); 

B=2*max(A(:,3)); 

gr=9.8; 
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[npts,d]=size(A); 

U=12.8; 

nfreq=500; 

T=5.5: 

beta=pi/3; 

F4=zeros(nfreq,1); 

H4=zeros(nfreq,1); 

for  nn=1:nfreq 
wo=2.5*nn/nfreq; 
k=wo''2/gr; 

we=wo-k*U*cos(beta); 

%  set  up  p  (position  matrix)  for  findfk 
for  sta=1:nsta 
clear  p 

q=1; 

for  r=1  :npts 
if  A(r,1)  ==  sta 
p(1,q)=A(r,2): 
p(2,q)=A(r,4); 
p(3,q)=A(r,3): 
q=q+1; 

x(sta)=A(r,2); 

end 

end 

%  get  sectional  FK,  Diffraction  Moments 
[f4(sta),h4(sta)]=findfk(p,  wo,  we,  k,T,  beta,B,sta); 
end 

%  Trapezoidal  Rule  Integration  of  Moments 
F4(nn)=0; 

H4(nn)=0: 
for  sta=1:nsta-1 
dx=x(sta)-x(sta+1); 

F4(nn)=F4(nn)+dx*(f4(sta)+f4(sta+1))/2; 

FI4(nn)=FI4(nn)+dx*(h4(sta)+h4(sta+1))/2; 

end 

F4(nn)=rho*gr*F4(nn); 

H4(nn)=rho*H4(nn); 

%  reset  for  next  frequency 
clear  f4 
clear  h4 

%  counter  for  status 
if  nn/10  ==  round(nn/10) 
fprintf('n=  %1.0f\n',nn) 
toe 
end 
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end 

toe 

1.a.  Subroutine:  ‘findfk.m’ 

function  [f4,  h4]  =  findfk(p,  w,  we,  k,  T,  beta,B,sta) 

%  Finds  sectional  Froude  Krylov  &  Diffraction  Excitation  Moments  for 
%  ForceFinder 

%  p  is  matrix  of  offsets,  a  is  incremental  waveheight,  psi  is  random  phase  angle, 
w,k  are  freq  &  wave  nr 

%KG  is  height  of  eg  above  baseline  T  is  still  water  draft 
%posit(1,:)  are  y-values,  posit(2,:)  are  z  values  in  the  section  plane 
%x  is  the  longitudinal  position  of  the  section 

%beta  is  the  angle  of  incidence  of  the  incoming  wave,  measured  in  radians  from 
the  stern 

%eta1  &  eta2  are  the  actual  waterlines  of  the  ship  after  translation  and 
application  of  the  wave 

%r  &  q  are  the  point  indices  corresponding  to  eta1  &  eta2 

%f2  &  f3  are  the  two-d  F-K  forces  for  the  section 

%posit  is  a  matrix  of  y,z  ordered  pairs 

posit(1,:)=p(2,:); 

posit(2,:)=p(3,:); 

x=p(1,1); 

[d,npts]=size(posit); 
posit(3,:)=zeros(1  ,npts); 
gr=9.8; 
c1=gr/w; 

c2=k*x*cos(beta); 

c3=k*posit(2,:)-T; 

%  This  block  of  code  calculates  the  2dimensional  diffraction  potential  for  the 

section 

npv=10; 

nph=55; 

npb=10; 

%  defines  domain  for  constant  panel  method 
[np,npanels,yvert,zvert,ybv,zbv,ycontrol,zcontrol,lgth,ny,nz]  = 
setpans(posit(1,:),posit(2,:)-T,npv,nph,npb,B,k,T); 

%  finds  G,  dG/dn  for  each  panel 

[g,  dgdn]  =  inflcoef(npanels,yvert,zvert,ycontrol,zcontrol,lgth); 

%  finds  A  &  b  in  discrete  BVP  A(phi)=b 
[Amatrix,b]  = 

findmatrix(np,nph,npv,npb,w,we,k,beta,x,ycontrol,zcontrol,ny,nz,g,dgdn); 

phid=Amatrix\b; 

%  determines  total  disturbance  potential 
for  m=1:np-1 

phil=c1*exp(i*(c2+k.*sin(beta)*posit(2,m)))*exp(c3(m)); 

PSI4(m)=phil+phid(np+nph-1+m); 
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phil=c1*exp(i*(c2-k.*sin(beta)*posit(2,m)))*exp(c3(m)); 

PSI4left(m)=phil+phid(np+nph-m); 

end 

%Translation  to  origin,  which  is  still  water  WL 
posit(2,:)=posit(2,:)-T; 

WL=0; 

f4=0; 

h4=0; 

%  q  is  number  of  points  below  waterline 
if  min(posit(2,:))>WL 
fprintf('station  is  above  water\n') 
return 
end 

if  max(posit(2,:))>WL 

q=1; 

while  posit(2,q)<WL 

q=q+1; 

end 

q=q-1; 

fprintf('%1 .2f  out  of  %1 .2f  points\n',q,npts-1 ) 
else 

q=npts-1; 

fprintf('station  is  below  water\n') 
end 
r=q: 

aa=zeros(3,npts); 

%find  forces  on  'the  right  half  of  the  hull 
for j=1:q 

%  compute  normals 

aa(:,j)=cross(posit(:,j+1)-posit(:,j),[0;0;-1]); 

amag=sqrt(aa(1  ,j)''2+aa(2,j)''2); 

n3(j)=aa(1,j)/amag; 

n2(j)=aa(2,j)/amag; 

mid=(posit(:,j)+posit(:,j+1))/2; 

n=[n2(j)  n3(j)  0]; 

n4v=cross(mid,n); 

n4(j)=n4v(3); 

z1=posit(2,j); 

z2=posit(2,j+1); 

dl=sqrt((posit(1  ,j)-posit(1  ,j+1  ))''2+(posit(2,j)-posit(2,j+1  ))''2); 

%  calculate  the  incremental  excitation  moment 

pt1  =exp(k*posit(1  ,j)*sin(beta)*i)*exp(k*z1  )*exp(-i*k.*x*cos(beta)); 

pt2=exp(k*posit(1,j+1)*sin(beta)*i)*exp(k*z2)*exp(-i*k.*x*cos(beta)); 
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f4=f4+n4(j)*(pt1+pt2)/2*dl; 

h4=h4-i*we*(n4(j))*dl*PSI4(j); 

end 

%find  forces  on  'the  left  half 

%first,  change  the  y-values  to  negative 
posit(1,:)=-posit(1,:); 

clear  n3 
clear  n2 
for  j=1:r 

aa(:J)=cross(posit(:,j+1)-posit(:,j),[0;0;-1]); 

amag=sqrt(aa(1  J)''2+aa(2J)''2); 

n3(j)=aa(1  j)/amag; 

n2(j)=aa(2J)/amag; 

mid=(posit(:,j)+posit(:J+1))/2; 

n=[n2(j)  n3(j)  0]; 

n4v=cross(mid,n); 

n4(j)=n4v(3); 

z1=posit(2J); 

z2=posit(2J+1); 

dl=sqrt((posit(1  J)-posit(1  J+1  ))''2+(posit(2J)-posit(2J+1  ))'^2); 

%  the  next  two  blocks  calculate  the  sum  of  discretized  waves  for 
%  exp(-ikxcosb)*exp(ikysinb)*e(kz),a=1 

pt1  =exp(k*posit(1  J)*sin(beta)*i)*exp(k*z1  )*exp(-i*k.*x*cos(beta)); 
pt2=exp(k*posit(1J+1)*sin(beta)*i)*exp(k*z2)*exp(-i*k.*x*cos(beta)); 
f4=f4+n4(j)*(pt1+pt2)/2*dl; 
h4=h4-i*we*(n4(j))*dl*PSI4left(j); 
end 

1.b.  Subroutine:  ‘setpans.m’ 

function  [np,npanels,yvert,zvert,ybv,zbv,ycontrol,zcontrol,lgth,nx,ny]  = 
setpans(yp,zp,npv,nph,npb,B,k,T) 


% 

%points  will  be  in  two  arrays,  yp  &  zp. 

%nph  is  #  of  horizontal  panels 

%npv  is  #  of  vertical  panels 

%np  is  #  points  in  arrays  yp  &  zp 

%hw  is  half  width  of  domain,  h  is  depth  of  domain 

%dl  is  length  of  top  hor.  panels,  dv  is  side  panels,  dw  is  bottom  hor  panels, 
if  max(zp)>0 
np=1; 
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d=0; 

while  zp(np)<0 
np=np+1; 
end 

yp=yp(1:np); 

zp=zp(1:np): 

dy=yp(np)-yp(np-1); 

dz=zp(np)-zp(np-1); 

ynew=yp(np)-zp(np)*dy/dz; 

yp(np)=ynew; 

zp(np)=0; 

end 

if  max(zp)<=0 
[d,np]=size(yp): 
end 

%  sets  the  domain  size:  hw  is  half-width  (10  beams) 

%  h  is  depth:  half  the  wavelength  or  three  times  the  draft, 

%  whichever  is  greater. 

%  for  the  ONR  hull,  B=18.8m  ->hw=188m,  nph=55,  dl=1.54m 


hw=10*B; 

h1=3/k; 

h2=3*T; 

h=max(h1,h2); 

dl=(hw-yp(np))/nph; 

dv=h/npv; 

dw=2*hw/npb; 

yvert(1)=-hw; 

zvert(1)=0; 

%  vertices  for  free  surface 
for  m=1:nph 
%  left  top  half 
yvert(m+1  )=yvert(m)+dl; 
zvert(m+1)=0; 

%  right  top  half 
trh=2*np+nph-1+m; 
yvert(trh)=yp(np)+dl*m; 
zvert(trh)=0; 
end 

%  vertices  for  hull 
counter=np; 
for  m=1:np 

yvert(m+nph)=-yp(counter); 

zvert(m+nph)=zp(counter); 

counter=counter-1; 
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mid=m+nph+np-1; 

yvert(mid)=yp(m); 

zvert(mid)=zp(m); 

end 

%vertices  for  sides 
for  m=1  :npv 

right=2*np+2*nph+m-1 ; 
left=right+npv+npb; 
yvert(right)=hw; 
zvert(right)=-m*dv; 
yvert(left)=-hw; 
zvert(left)=-h+m*dv; 
end 

%  vertices  for  bottom 
for  m=1:npb 

bot=m+2*np+2*nph+npv-1 ; 
yvert(bot)=hw-dw*m; 
zvert(bot)=-h; 
end 

endpt=2*(np+npv+nph)+npb-1 ; 
yvert(endpt)=yvert(1 ); 
zvert(endpt)=zvert(1 ); 
npanels=endpt-1; 

%  control  points 
for  k=2  :endpt; 

ycontrol(k-l)  =  0.5*(yvert(k-1)+yvert(k)); 
zcontrol(k-l)  =  0.5*(zvert(k-1)+zvert(k)); 
lgth(k-1  )=sqrt((yvert(k)-yvert(k-1  ))'^2  ... 

+(zvert(k)-zvert(k-1  ))''2); 
nx(k-1  )=  -(zvert(k)-zvert(k-l  ))/lgth(k-1 ); 
ny(k-1)  =  (yvert(k)-yvert(k-1))/lgth(k-1); 
end 

nbv=endpt-(2*np-3); 

ybv(1  :nph+1  )=yvert(1  :nph+1 ); 

zbv(1  :nph+1  )=zvert(1  :nph+1 ); 

ybv(nph+2:nbv)=yvert(2*np+nph-1;endpt); 

zbv(nph+2:nbv)=zvert(2*np+nph-1  :endpt); 


1.C.  Subroutine:  ‘findmatrix.m’ 

function  [Amatrix,b]  = 

findmatrix(np,nph,npv,npb,w,we,k,beta,x,yc,zc,ny,nz,g,dgdn) 
%  X  is  longitudinal  position  of  station 
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%  initialize  variables 
gr=9.81; 
c1=i*gr/w; 
c2=k*x*cos(beta); 

npanels=2*nph+2*np+2*npv-2+npb; 

b=zeros(npanels,1); 

Amatrix=zeros(npanels,npanels); 

for  m=1:npanels 
for  n=1:nph 

%  top  left  surface;  w''2phi-gdphi/dn=0; 

Amatrix(m,n)=dgdn(m,n)-w''2/gr*g(m,n); 

%top  right  surface;  same  boundary  condition 
nn=n+2*np+nph-2; 

Amatrix(m,nn)=dgdn(m,nn)-w''2/gr*g(m,nn); 

end 

for  n=nph+1:2*np+nph-2 

%ship's  hull;  dphi/dn=-dphl/dn 
Amatrix(m,n)=dgdn(m,n); 

phil=c1*exp(i*(c2-k*sin(beta)*yc(n-nph))+k*zc(n-nph)); 

vv=-phiri*k*sin(beta); 

ww=k*phil; 

vdotn=vv*ny(n-nph)+ww*nz(n-nph); 

b(m)=b(m)-g(m,n)*vdotn; 

end 

for  n=1;npv 

%right  side  dphi/dn=-ikphi 
nn=n+2*nph+2*np-2; 

Amatrix(m,nn)=dgdn(m,nn)+i*k*g(m,nn); 

%left  side  same  boundary  condition 
nn=n+2*nph+2*np-2+npv+npb; 
Amatrix(m,nn)=dgdn(m,nn)+i*k*g(m,nn); 
end 

for  n=1:npb 

%bottom  dphi/dn=0 
nn=n+2*nph+2*np-2+npv; 

Amatrix(m,nn)=dgdn(m,nn); 

end 

Amatrix(m,m)=Amatrix(m,m)+pi; 

end 

2.  Main  Program:  LinearResponse  This  program  only  needs  the 
raofunctions.mat  file  to  run. 

%  Calculates  Time  simulations  for  heave,  pitch,  and  roll 


59 


%  Needs  raofuncations.mat  to  run 
%  raofunctions  has  the  following  row  vectors; 

%  freqs  (wave  frequencies  of  interest) 

%  efreq  (encounter  frequencies  at  each  wave  frequency) 
%  A  (wave  amplitudes  at  each  wave  frequency) 

%  RRAO  (roll  RAO) 

%  Rph  (roll  phase  in  degrees) 

%  HRAO  (heave  RAO) 

%  Hph  (heave  phase  in  degrees) 

%  PRAO  (pitch  RAO) 

%  Pph  (pitch  phase  in  degrees) 

%  save  output  as  linresplOO 

clear 
close  all 

load  raofunctions 

tic 

[b,a]=size(Hph); 

Hph=Hph*pi/180; 

Pph=Pph*pi/180; 

Rph=Pph*pi/180; 

k=freqs.''2/9.8; 
we= efreq; 

%  allow  for  more  than  one  run  for  statistics 

runs=1; 

tinc=100; 

dt=1/tinc; 

%  previous  two  lines  set  dt  to  1/tinc  seconds 
npts=3600*tinc; 
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t=(1:npts)*dt; 
heave=zeros(npts,1 ); 
pitch=zeros(npts,1); 
roll=zeros(npts,1); 
wh=zeros(npts,1); 
for  n=1  :runs 

%  set  new  random  phase  angles  for  discretized  waves 
rand('state',sum(1 00*clock)); 
psi=2*pi*rand(length(A),1 )'; 
for  m=1:npts 

wh(m)=sum(A.*cos(efreq*t(m)+psi)); 
heave(m)=sum(A.*HRAO.*cos(abs(efreq)*t(m)+psi+Hph)); 
pitch(m)=sum(k.*A.*PRAO.*cos(abs(efreq)*t(m)+psi+Pph)); 
roll(m)=sum(k.*A.*RRAO.*cos(abs(efreq)*t(m)+psi+Rph)); 
if  m/tinc==round(m/tinc) 
time=m/tinc 
toe 
end 
end 

draft=5.5+wh-heave; 

maxdraft(n)=max(draft); 

mindraft(n)=min(draft); 

n 

toe 

end 


3.  Main  Program:  rollintegration: 

%  Integrates  exciting  forces  over  time  to  produce  a  time  sim  of  roll  angle 
%  Needs  RightingMoment.m  to  run 
%  Needs  the  following  files  to  run: 

%  Rmomentdata.mat; 

%  Contains  data  to  conduct  the  lookup  table  function  for  righting  moment 
%  forcedataSOO.mat:  Output  of  ForceFinder,  saved  as  a  .mat  file 
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%  linresplOO.mat:  Output  of  LinearResponse,  saved  as  .mat  file 


clear  all 
close  all 

global  RM 
global  rolls 
global  drafts 
global  pitchs 
tic 

load  Rmomentdata 
load  forcedataSOO 
load  linresplOO 

%  sets  the  time  simulation  to  npts*dt  seconds 

npts=360000; 

dt=.01; 

%  Sets  the  moment  of  inertia,  added  inertia,  damping 
%  and  linear  restoring  coefficients 
mass=8402180; 
lactual=mass*(.4*B)''2; 

A44=.3*lactual; 

l=lactual+A44; 

b=.075; 

c44=1 .476*mass*gr; 

B44=b*sqrt(c44*l); 

%  initialize  the  catastrophic  roll  event  counter 
count=0; 

theta=zeros(npts,1); 
thdot=zeros(npts,1); 
thddot=zeros(npts,1 ); 
theta2=zeros(npts,1 ); 
thdot2=zeros(npts,1 ); 
thddot2=zeros(npts,1 ); 

M=zeros(npts,1); 

C=zeros(npts,1); 

C2=zeros(npts,1); 

%  provide  the  option  to  trim  by  the  stern 
pitch=pitch; 

%  set  IC's  to  match  linresplOO 

theta(1)=roll(1); 

thdot(1)=(roll(2)-roll(1))/dt; 

theta2(1)=theta(1); 

thdot2(1)=thdot(1); 
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%Set  magnitude  and  phase  of  excitation  forces 

exmoment=abs((F4+H4)'.*A); 

mphase=angle(F4+H4)'; 

for  m=1:npts-1 

M(m)=sum(exmoment.*cos(we*t(m)+psi+mphase)); 

%  protect  for  high  pitch  values 
if  pitch(m)>max(pitchs) 
pitch(m)=max(pitchs); 
end 

%  check  for  catastrophic  roll  event 
if  theta(m)>pi/2 
theta(m)=theta2(m); 
thdot(m)=thdot2(m); 
count=count+1 
end 

if  theta(m)<-pi/2 
theta(m)=theta2(m); 
thdot(m)=thdot2(m); 
count=count+1 
end 

%  Find  non-linear  restoring  moment 
C(m)=9.8*RightingMoment(draft(m),pitch(m),abs(theta(m))); 
if  theta(m)<0 
C(m)=-C(m): 
end 

%  Find  linear  restoring  moment 
C2(m)=c44*theta2(m); 

%  Forward  Euler  integration 
thddot(m)=(M(m)-C(m)-B44*thdot(m))/l; 
thdot(m+1)=thdot(m)+thddot(m)*dt; 
theta(m+1)=theta(m)+thdot(m)*dt; 

thddot2(m)=(M(m)-C2(m)-B44*thdot2(m))/l; 

thdot2(m+1)=thdot2(m)+thddot2(m)*dt; 

theta2(m+1)=theta2(m)+thdot2(m)*dt; 

end 

t=t/60; 

figure(1) 

plot(t(1  :m),theta2(1  :m)*1 80/pi, t(1  :m),roll(1  :m)*1 80/pi) 
figure(2) 

plot(t(1  :m),theta(1  :m)*1 80/pi, t(1  :m),theta2(1  :m)*1 80/pi) 
toe 
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%  function  to  determine  non-linear  roll  moment 
%  variable  RM  is  lookup  table  for  righting  moments. 
%  RM(:,p,r)  is  moment  vs.  draft 
%  RM(d,:,r)  is  moment  vs.  pitch 
%  RM(d,p,:)  is  moment  vs.  roll 
%  From  an  excel  spreadsheet: 

% 

%  roll:  0  .05  ... 

%  pitch 

%  -.15  RM(d,-.15,0)  RM(dr.15,.05) 

%  -.1  RM(d,-1,0)  RM(d,-.1,.05) 

% 


3.a  subroutine:  ‘RightingMoment’ 

function  moment  =  RightingMoment(heave,  pangle,  rangle) 

global  RM 
global  drafts 
global  pitchs 
global  rolls 

%  protects  from  high  pitch/heave 
if  heave<2 
heave=2 
end 

if  pangle<-0.1337 
pangle—. 1337 
end 

rangle=abs(rangle); 

%  find  heave  index 
d=1; 

while  drafts(d)<=heave 
d=d+1; 
end 

%  find  pitch  index 
P=1; 

while  pitchs(p)>=pangle 
P=P+1; 
end 

%  find  roll  index 
r=1; 

while  rolls(r)<=rangle 
r=r+1; 
end 


%  3-D  interpolation 


64 


HMRA=RM(d,p,r)-(rolls(r)-rangle)/(rolls(r)-rolls(r-1))*(RM(d,p,r)-RM(d,p,r-1)); 
LMRA=RM(d,p-1  ,r)-(rolls(r)-rangle)/(rolls(r)-rolls(r-1  ))*(RM(d,p-1  ,r)-RM(d,p-1  ,r- 
1)); 

HRA=HMRA-(pitchs(p)-pangle)/(pitchs(p)-pitchs(p-1))*(HMRA-LMRA); 
HMRA=RM(d-1  ,p,r)-(rolls(r)-rangle)/(rolls(r)-rolls(r-1  ))*(RM(d-1  ,p,r)-RM(d-1  ,p,r- 
1)); 

LMRA=RM(d-1  ,p-1  ,r)-(rolls(r)-rangle)/(rolls(r)-rolls(r-1  ))*(RM(d-1  ,p-1  ,r)-RM(d-1  ,p- 
1,r-1)); 

LRA=HMRA-(pitchs(p)-pangle)/(pitchs(p)-pitchs(p-1))*(HMRA-LMRA); 

moment=HRA-(drafts(d)-heave)/(drafts(d)-drafts(d-1))*(HRA-LRA); 

return 
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Appendix  B:  List  of  Symbois 


Variable 

Meaning 

A 

Wave  Amplitude 

Ajk 

Added  Mass 

B 

Beam 

Bjk 

Damping  Coefficient 

Cjk 

Restoring  Coefficient 

f 

Section  Froude-Krylov 
Moment 

Fj 

Force  In  The  j  Direction 

g 

Acceleration  Of  Gravity 

G 

Green’s  Function 

h 

Sectional  Diffraction 
Moment 

H 

Transfer  Function 

H4 

Diffraction  Moment 

H1/3 

Significant  Wave  Fleight 

I44 

Roll  Moment  Of  Inertia 

k 

Wave  Number 

k44 

Roll  Gyradius 

i 

Length 

Mjk 

Generalized  Mass  Term 

Hi 

Normal  Vector 

Component  In  J 

Direction 

p 

Pressure 

RAO 

Response  Amplitude 
Operator 

S 

Wave  Spectral  Density 

t 

Time 

u 

Speed  In  The  X- 
Direction 

Variable 

Meaning 

V 

Speed  In  The  Y- 
Direction 

V 

Velocity  Vector 

V 

Volume 

w 

Speed  In  The  Z- 
Direction 

X 

Longitudinal  Direction 

X44 

Non-dimensional  Added 
mass  coefficient 

y 

Transverse  Direction 

z 

Vertical  Direction 

a 

Phase  Angle  For 

Transfer  Function  or 
Moment 

3 

Wave  Angle  Of 

Incidence 

P44 

Non-Dimensional 
Damping  Term 

A 

Displacement 

C 

Wave  Elevation 

n 

Response  Value 

0 

Roll  Angle 

A 

Tuning  Coefficient  Or 
Wavelength 

p 

Density 

0 

Velocity  Potential 

0D 

Diffraction  Potential 

01 

Incident  Potential 

0 

Wave  Phase  Angle 

O) 

Circular  Frequency 

We 

Encounter  Frequency 

Wm 

Modal  Frequency 

Wo 

Wave  Frequency 
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